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ABSTRACT 

Newtonian simulations are routinely used to examine the matter dynamics on non¬ 
linear scales. However, even on these scales, Newtonian gravity is not a complete 
description of gravitational effects. A post-Friedmann approach shows that the leading 
order correction to Newtonian theory is a vector potential in the metric. This vector 
potential can be calculated from N-body simulations, requiring a method for extracting 
the velocity field. Here, we present the full details of our calculation of the post- 
Friedmann vector potential, using the Delauney Tesselation Field Estimator (DTFE) 
code. We include a detailed examination of the robustness of our numerical result, 
including the effects of box size and mass resolution on the extracted fields. We present 
the power spectrum of the vector potential and find that the power spectrum of the 
vector potential is ~ 10 5 times smaller than the power spectrum of the fully non¬ 
linear scalar gravitational potential at redshift zero. Comparing our numerical results 
to perturbative estimates, we find that the fully non-linear result can be more than an 
order of magnitude larger than the perturbative estimate on small scales. We extend 
the analysis of the vector potential to multiple redshifts, showing that this ratio persists 
over a range of scales and redshifts. We also comment on the implications of our results 
for the validity and interpretation of Newtonian simulations. 

Key words: gravitation - cosmology: theory - cosmology: large-scale structure of 
the universe. 


On the largest scales in cosmology, theoretical calcula¬ 
tions can be carried out using standard cosmological pertur¬ 
bation theory. These calculations fully encompass General 
Relativity (GR) but are limited to scales where the pertur¬ 
bations, in particular the density perturbation, are small. On 
smaller scales, where the focus is on non-linear structure for¬ 
mation, Newtonian N-body simulations are used. These sim¬ 
ulations do not require that the density contrast is small, but 
they suffer from the limitations of being Newtonian rather 
than GR simulations. There is an entire field in cosmology 
dedicated to developing, running and analysing these New¬ 
tonian N-body simulations. There has been sporadic interest 
in understanding the use of Newtonian theory in cosmology 


(Tomita||1991 

Shibata & Asada||1995 Takada & Futamase 

1997 

|Matarrese & Terranova 1996 Carbone & Matarrese 

2005 

|Hwang et al.||2008; Hwang & Noh 2013 Milillo et al. 


2015 

Flender & Schwarz120121 Haugg et al. 2012 Kopp et al.| 

2014 

1, as well as examining the relativistic interpretation of 

the simulations (Chisari & Zaldarriaga|2011 Green & Wald 

2012 

Bruni et al. 2014 Adamek et al. 2013). These stud- 


ies have predominantly focussed on whether the dynamics of 
density contrast and scalar potential accurately match those 


of GR. 

In this paper, we are mostly interested in another im¬ 
portant limitation of Newtonian simulations: Even if the 
matter dynamics are being computed correctly, there are 
cosmological quantities of interest on non-linear scales that 
have no counterpart in Newtonian theory. Examples of these 
quantities include the difference between the two scalar po¬ 
tentials, gravitational waves and the vector potential in the 
metric, all of which must exist on non-linear scales in a GR 
universe. 
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These extra quantities would naively be expected to be 
small if the Newtonian simulations are a good approxima- 
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tion to a GR universe. However, explicitly calculating these 
quantities has several advantages: To start with, it would be 
good to have a quantitative check of whether these quanti¬ 
ties are small, and indeed how small they are. In particular, 
as we enter the era of precision cosmology, we need to check 
that these quantities will not affect the observables at the 
percent-level. Furthermore, checking that these quantities 
are negligible provides a quantitative check on the Newto¬ 
nian approximation in a ACDM cosmology. 

We will be working with the post-Friedmann formalism 
| Milillo et al.|2015 Milillo|2010 I . This generalises to cosmol¬ 
ogy the weak-field (post-Minkowski) approximation, with 


a post-Newtonian style expansion 1 

Chandrasekhar 1965 

Weinberg 1972 Poisson & Will|2014 

) in inverse powers of 


the speed of light c of the perturbative quantities. These 
expansions need to be performed differently in cosmology 
compared to in the Solar System due to the different sit¬ 
uations and aims in the two cases. For example, the time¬ 
time and space-space components of the metric need to be 
treated at the same order in cosmology in order for the re¬ 
sulting equations to be a consistent solution of the Einstein 
equations. 

The post-Friedmann formalism, when linearised, cor¬ 
rectly reproduces conventional linear perturbation theory 
and can thus describe structure formation on the largest 
scales. More importantly, the leading order equations in the 
1/c expansion can be examined and are expected to yield the 
non-linear Newtonian equations. Note that in this “Newto¬ 
nian” regime, the density contrast has not been assumed to 
be small. The equations in this regime will be shown in sec¬ 
tion [T] essentially comprising the Newtonian equations, as 
expected, plus an additional equation. This additional equa¬ 
tion shows how the vector potential in the metric, the lowest 
order beyond-Newtonian quantity, is generated by the mat¬ 
ter dynamics. This vector potential is the beyond-Newtonian 
quantity that we will examine in this paper, it is the cosmo¬ 
logical manifestation of the ubiquitous relativistic effect of 
frame dragging. This effect has been measured in the Solar 
System by Gravity Probe B (Everitt et al|201ip . 

In this paper, we present a calculation of this vector 
potential based on extracting the density and velocity fields 
from N-body simulations. We expand on the results of |Bruni| 


et al. (20141, which was the first calculation of an intrin¬ 


sically relativistic quantity on fully non-linear scales from 
large scale cosmological matter fields, rather than from in¬ 
dividual astrophysical occurences. The main focus in this 
paper is to present the method used to extract this vector 
potential from N-body simulations. In particular, we exam¬ 
ine the robustness of the numerical extraction of the vector 
potential and present the tests we carried out to examine the 
numerical effects of simulation parameters on the extraction, 
which were not presented in Bruni et al. (2014). 

The main physical results of this paper are figures [3] 
and [8] showing the power spectrum of the vector potential 
at redshift zero and its evolution with time respectively. Ad¬ 
ditionally, we have presented the ratio of the vector potential 
power spectrum to that of the scalar potential in figures [5] 
and [9] The results on the magnitude and evolution of the 
power spectrum of the vector potential in this paper were 
used in |Thomas et al. (20141 to examine the possible weak- 
lensing consequences of the vector potential. 

This paper is laid out as follows. In section[l] we present 


the pertinent details of the post-Friedmann formalism and 
show the equation governing the vector potential. We will 
also present our definitions and notation regarding vector 
power spectra. In section [2] we explain how the relevant 
fields were extracted from N-body simulations and examine 
the robustness of this extraction. In section [3] we show the 
power spectrum of the vector gravitional potential and its 
time evolution, as well as comparing it to the closest an¬ 
alytical results in the literature. We conclude in section [4] 
Appendix [X] contains some details about vector power spec¬ 
tra and in Appendix [B] we show results from some of the 
numerical tests that were carried out. Additional plots are 
available as online supplementary material, divided amongst 
three files: Resolution_and_BoxSize_Dependence.pdf (here¬ 
after RB), GridSize_and_Binning_Dependence.pdf (hereafter 
GB) and ConsistencyChecks.pdf (hereafter CC). 


1 POST-FRIEDMANN FORMALISM 


The post-Friedmann approach is developed in |Milillo et al.| 
(20151 and Milillo (20101, see there for the full details. This 


approach considers a dust (pressure-less matter) cosmology 
with a cosmological constant. The perturbed FLRW metric, 
in Poisson gauge, is expanded up to order c -5 , keeping the 
goo and gij scalar potentials at the same order: 


goo = - 

got = 

gij ~ 


2Un 


1 




aBf 


aBj 


(1) 


1 + -Jf + ( 2V N + 4 Vp) 


x.. + 


The goo and g,:j scalar potentials have been split into the 
Newtonian ( Un , Vn) and post-Friedmann (Up, Vp) com¬ 
ponents. Similarly, the vector potential has been split up 
into Bf and Bf. Since this metric is in the Poisson gauge, 
the three-vectors Bf and Bf are divergenceless, Bf = 0 
and Bf, = 0. In addition, hij is transverse and tracefree, 
h\ = hfj = 0. Note that at this order, hij is not dynam¬ 
ical, so it does not represent gravitational waves. From a 
post-Friedmann viewpoint, there are two different levels of 
perturbations in the theory, corresponding to terms of order 
c ~ 2 and c -3 , or of order c -4 and c -5 respectively. Defining 
“resummed” variables, such as $ = 2Un+c~ 2 (2Uf — 4 Up), 
then calculating the Einstein equations and linearising them, 
reproduces linear GR perturbation theory in Poisson gauge. 
Thus, this approach is capable of describing structure for¬ 
mation on the largest scales. 

For smaller scales, in a dust cosmology, we are interested 
in the weak field, slow motion, sub-horizon, quasi-static and 
negligible pressure regime. This is simply derived by retain¬ 
ing only the leading order terms in the c _1 expansion and 
upon doing so we recover Newtonian cosmology, albeit with 
a couple of subtleties. The first is that the space-time met¬ 
ric is a well-defined approximate solution of the Einstein 
equations. The second is that we have an additional equa¬ 
tion, which is a constraint equation for the vector gravita¬ 
tional potential Bf. The full system of equations obtained 
from the Einstein and hydrodynamic equations ( |Mil illo et al. 
2015), given the evolution of the background a(t), is as fol- 
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lows. 


dvi a 1 _ _ 

—77—I- Vi = — UN,i 

at a a 

1 V7 2 T/ 47rG - 

V Uv =- ^pbS 

c z a z c z 

-^-V 2 (V N - U N ) = 0 

c z a z 

2d TT 2- 1 

+ a^ N,i - 2a? 


V~Bi 


8-irGpb 


( 2 ) 

(3) 

(4) 

(5) 

(1 + S) v(6) 


As expected, we have the Newtonian continuity and Euler 
equations from the hydrodynamic equations as well as Pois¬ 
sons equation from the Einstein equations. Note that the 
time derivative here is the convective derivative, dA/dt = 
dA/dt + v l A^/a, for any quantity A. The Einstein equations 
yield two additional equations: The first is an equation forc¬ 
ing the scalar potentials Vi v and Un to be equal, consistent 
with there being only one scalar potential in Newtonian the¬ 
ory. Note that some approaches consider the potentials to 
be a priori equal at leading order whereas here we assumed 
the full generality of GR and the equality of the potentials 
arose naturally on taking the Newtonian regime. The second 
additional equation relates the leading order vector gravita¬ 
tional potential, , to the momentum of the matter. Thus, 
even in the regime where the matter dynamics are correctly 
described by Newtonian theory, the frame-dragging poten¬ 
tial Bf 1 should not be set to zero; this would correspond 
to putting an extra constraint on the Newtonian dynamics. 
We note that there is a similar equation in several other for¬ 
malisms in the literature (Takada fe Futamase|1997 Hwang 


& No h|2013 Green & Wald 2012 I. We can see from equation 
^ that the potential B t is sourced by the vector part of 
the energy current pv. This is made apparent by taking the 
curl of this equation, which gives 


V x V 2 B n = - (l6?vGp b a 2 ) V x [(1 + 5)v\, 


(7) 


where the source term on the right hand side splits up into 
three terms: the vorticity Vxf and then two further terms, 


V x [(1 + <5)F] =Vxf + <SV x v + VS x v. 


( 8 ) 


It is equation 0 that will be used for the rest of the 
paper. Since the matter dynamics are not affected at this 
order, i.e. they are described by the standard Newtonian 
equations ( 2)4 1, the density and velocity fields sourcing the 
vector potential are Newtonian and can be extracted from 
N-body simulations. Using the definitions of vector power 
spectra in Appendix [XJ the power spectrum of the vector 
potential is given by 


P Mk)= ( 16 ^r a2 ) 


(9) 


with 


Psv = Pvxv{k) + Psvxv{k) + P(-vs)xv{k) (10) 

+-P(V<Sxv)(VxC)(&) + PcV5xv)(6VXv){k) + P(SS7 X v)(V X v) (k) ■ 

Unless stated otherwise, all plots of the gravitational po¬ 
tentials show the dimensionless power spectrum A(fc), see 
Appendix [A] for conventions. 


Table 1. Parameters for the simulations. 


Box size 
(/i -1 Mpc) 

Particle 

number 

Mass resolution 
(10 s Mo) 

Number of 
Realisations 

Softening 
(/i -1 kpc) 

80 

512 3 

3.97 

8 

6.25 

80 

512 3 

3.97 

1 

4.0 

140 

768 3 

6.31 

8 

6.25 

140 

560 3 

16.3 

8 

6.25 

160 

1024 3 

3.97 

3 

6.25 

160 

880 3 

6.26 

3 

6.25 

160 

640 3 

16.3 

8 

6.25 

160 

640 3 

16.3 

1 

5.0 

160 

320 3 

130 

8 

15.0 

200 

1024 3 

7.76 

2 

6.25 

240 

960 3 

16.3 

3 

6.25 

240 

480 3 

130 

8 

15.0 

320 

640 3 

130 

8 

15.0 


Table 2. Redshifts used to probe time evolution of quantities. 


Scale Factor 

Redshift 

Colour on time evolution plots 

0.33 

2.0 

black 

0.4 

1.5 

red 

0.5 

1.0 

magenta 

0.6 

0.67 

yellow 

0.7 

0.43 

green 

0.8 

0.25 

cyan 

0.9 

0.11 

blue 

1.0 

0.0 

brown 


2 SIMULATIONS 


Our simulations have all been run using the publicly avail¬ 


able N-body code Gadget2 (Springel 20051. Many simula¬ 


tions have been run in order to quantify the effects of box 
size and mass resolution on the quantities that we are ex¬ 
tracting, see table [l] for a full list of the simulations. All of 
the simulations were run with dark matter particles only, as 
the equation for the vector potential is derived for a pres¬ 
sureless matter and cosmological constant cosmology. To al¬ 


low comparison to previous studies of vorticity (Pueblas & 
|Scocci marro 2 009[ ), the simulations were run with the cos¬ 
mological parameters fl m = 0.27, SIa = 0.73, = 0.046, 

h = 0.72, r = 0.088, as = 0.9 and n s = 1. All of the simula¬ 
tions started at redshift 50 and had their initial conditions 
created using 2LPTic (Crocce et ah|2006). Our final result 
for the vector potential is taken from the three 160h~ 1 Mpc 
simulations with 1024 3 particles, these will be referred to as 
the high-resolution (HR) simulations. 


2.1 Tesselation 

To extract the necessary fields from the simulations, the De- 
launey Tesselation Field Estimator (DTFE) code was used 
(Cautun & van de Weygaert|201 1). Standard methods of ex¬ 
tracting fields from N-body simulations, such as Cloud-In- 
Cells (CIC) ( jHockney fc Eastwood| 1 981| ) work well for the 
density field, as the particles, by definition, sample the den¬ 
sity field well. However, these methods have several short¬ 
comings when applied to the extraction of velocity fields: 
One is that the field is only sampled where there are par- 
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tides, so in a low density region the velocity field is arti¬ 
ficially set to zero. In addition, the extracted field will be 
a mass-weighted, rather than volume-weighted field. A con¬ 
sequence of these shortcomings is that, as the grid size is 
increased, the velocity field will not converge. In fact, it will 
become zero in an increasing proportion of the grid cells as 
the grid size increases. Several authors have looked at us¬ 


ing the Delauney tesselation (Schaap & van de Weygaert 


2000[ [van de Weygaert fe Schaap||2009| |Bernardeau fc van 

de Weygaert||1996|) for astrophysical applications including 

the examination of velocity fields. See also |Pueblas fe Scoc-| 
(2009) for comparisons of extracting velocity fields 


with tesselations rather than more standard methods. The 
DTFE code constructs the Delauney tesselation of the set of 
particles, consisting of tetrahedra whose nodes are located 
at the particles positions. The tetrahedra are constructed 
such that the circumsphere of each tetrahedron does not 
contain any of the particles except for the particles located 
at the nodes of the tetrahedron in question. This makes the 
tesselation unique. The particles’ velocities are then linearly 
interpolated across each tetrahedron, yielding a value for the 
smoothed velocity field and its gradients at every point in 
the simulation volume. A regular N| rid grid is laid down and 
the code samples Ai aamp i ea points at random in each grid cell 
and averages the field over these points, giving a value for the 
smoothed field in each grid cell. Once the fields are obtained 
on the regular grid, the power spectra are calculated using 
the standard process of averaging the modulus-squared of 
the Fourier coefficients over a given range of k. For the anal¬ 
yses here, we used Al gr id/4 bins, although varying this value 
does not affect the results (see Appendix B71. 


2.2 Convergence and Tests 

It is important to ensure that our numerical result for the 
vector potential is robust and independent of the simulation 
parameters. In this subsection we will present the results of 
our examination into the effects of different simulation pa¬ 
rameters on the extracted vector power spectrum. Since the 
velocity and density fields both contribute to the source for 
the vector potential, we will examine the density, vorticity 
and velocity divergence spectra too: We will examine their 
behaviour individually, compare them to other studies and 
methods of extraction and also consider the consistency of 
the extracted fields through the relations 

k 2 P s {k ) = P vs (k) 

k 2 Pv{k) = Pv-v(k) + Pvxv{k). (11) 

The box size and mass resolution of the simulation are 
the two main parameters whose effect on the extracted fields 
needs to be examined. In addition, we have examined the 
effect of varying the grid size and Ai a am P iea, which are both 
internal DTFE parameters. The parameters of the differ¬ 
ent simulations used are in table [TJ We chose the soften¬ 
ing lengths of the N-body simulations to be consistent with 
Puebla s" fe Scoccimarro| (| 2009) in order to recreate their 
study of the velocity divergence and vorticity, however vary¬ 
ing the softening length did not influence the results, see 
Appendix |B5| 

Although we did run some simulations with a box size 
below 140/i _1 Mpc, we have not included these in the analy¬ 


sis here as smaller box sizes have systematically less power. 
See Appendix |B6| for the results from these simulations and 
how they compare to the larger box sizes. For further results 
regarding the effects of a small simulation box on cosmologi¬ 
cal quantities, see (Bagla & Prasad 2006 Bagla fe Ray|2005 


Gelb & Bertschinger|]19941. 


2.2.1 A note on error bars 

Since we have only three realisations of our HR simulations, 
we cannot compute meaningful error bars. Thus, we have 
not included any error bars on the majority of our plots. 
Instead, in figures [l] [2] [5] and [6j we have plotted the results 
from the three individual realisations, in order to illustrate 
by how much the results vary. Unless stated otherwise, the 
results shown in the other plots show the average over the 
realisations. We explicitly examine the variation amongst re¬ 
alisations in Appendix [B4] for several quantities, notably the 
vorticity and vector potential. In particular, we note there 
that when considering the vector potential, cosmic variance 
on the largest scales affects smaller scales, as explained by 
a perturbative analysis ( |Lu et al.|2008| |Hui-Ching Lu et al.| 
20091. See B4 for more discussion of this. We also note there 


that the variation of the vorticity amongst realisations seems 
to be larger than the variation of the density, although there 
seems to be no discussion of this in the literature. 


2.2.2 Mass Resolution 

We have examined the dependence of the density, veloc¬ 
ity divergence, vorticity and vector potential on the mass 
resolution of the simulations. For the density and velocity 
divergence there is evidence for a mild dependence on mass 
resolution for both of these fields on smaller scales. This 
is likely to be due to the DTFE window function, which 
cannot be compensated for, rather than a mass-resolution 
dependence of the field itself. There is no evidence of any 
mass-resolution dependence of these fields on larger scales. 
The variation of the density and velocity divergence with 
mass resolution can be seen in figures 1 and 2 of file RB. 
The effect of the small-scale mass-resolution dependence is 
negligible for our HR simulations, as seen when compar¬ 
ing to alternative methods of calculating the density power 
spectrum. 

The dependence of the vorticity power spectrum with 
mass resolution is shown in figure [Bl] The power spectrum 
shows spurious additional power when the mass resolution 
is insufficient. However, once the resolution is sufficient, of 
order 10 9 Mq, there is no evidence for any systematic de¬ 
pendence on mass resolution. This dependence on mass res¬ 
olution, followed by convergence around ~ 10 9 Mq, matches 
previous findings, notably those of|Pueblas & Scoccimarro| 
( 2009 1. 

In figure |B2] we show the dependence of the vector po¬ 
tential on mass resolution. There is a clear dependence of 
the vector potential on mass resolution, similar to that seen 
for the vorticity. However, there are several differences. In 
particular, the mass-resolution dependence seems to be less 
important for smaller scales, where there is a greater depen¬ 
dence on box size (see later). In addition, the dependence 
on mass-resolution is still apparent around 10 9 Mq. How¬ 
ever, once there mass resolution has improved to around 
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6 x 10 8 Mq, there is no evidence of a mass resolution depen¬ 
dence of the vector potential. 

To show this further, figure |B8| shows the higher res¬ 
olution simulations in more detail, complete with the indi¬ 
vidual realisations of the HR simulations. The y-axis here 
is k 2 Pg{k) in order to show the variance more clearly over 
the range of scales being considered. The cyan line shows 
the simulation with the worst resolution (16.3 x 10 8 Mq) of 
those in this plot and indeed this simulation shows a system¬ 
atic deviation on the largest scales. The better resolution 
simulations show better convergence, with the 140h _1 Mpc 
simulation with 768 3 particles being consistent with the HR 
simulations for essentially the entire range under considera¬ 
tion. This convergence is examined further in Appendix |B4| 

2.2.3 Box Size 

We have also considered the effect of varying the box size 
on the extracted power spectra. As expected, there is no 
evidence for any systematic dependence of the density, vor- 
ticity and velocity divergence power spectra on the box size 
of the simulations. This can be seen in figures 3, 4 and 5 of 
file RB. Note that, for sufficiently small boxes, a systematic 
deviation can arise, see Appendix |B6| 

Figure [B3| shows the box size dependence of the vector 
potential. As mentioned above, the vector potential does 
show some dependence on box size. The vector potential 
shows signs of a dependence on the box size on scales below 
l/i -1 Mpc, however this is difficult to entangle from the ef¬ 
fects of mass resolution and the window function. For box 
sizes below 200/iT 1 Mpc, there is no systematic dependence 
of the vector potential power spectrum with box size. 

In Appendix |B4| we examine the variation between realisa¬ 
tions for the vector potential, and relate it to the behaviour 
that might be expected from perturbative arguments. In 
particular, figure |B8| shows how the variation between re¬ 
alisations is larger than the effects of box size and mass 
resolution for simulations with box sizes below 200/i _1 Mpc 
and mass resolution of at least 6 x 10 8 Mq. Thus, we expect 
numerical effects from the simulation parameters to be a 
sub-dominant source of error as long as the parameters are 
within this range. 

2.2.4 Consistency Checks 

There are a few consistency checks that can be performed on 
the different fields that we are interested in. The quantities 
that are used for the vector potential include the density 
field and its gradients as well as the velocity field and its 
gradients. There are two relations between these fields and 
their derivatives, 

k 2 P s {k ) = P vs (k) (12) 

k 2 Pv{k) = Pv-v(k) + Pvxv{k). (13) 

We have extracted the quantities on the left and right sides 
of these relations from our HR simulations and compared 
them, see figure 1 in file CC for the ratio P\rs(k)/k 2 Ps{k) 
and figure 2 in file CC for the ratio k 2 P$(k )/(Pv ff(fc) + 
P\7xv{k)). In both cases, two curves are plotted, correspond¬ 
ing to two different methods of calculating the ratio. The 
blue line shows the ratio exactly as suggested above, with 


the factor of k in equation |12| taken to be the value defining 
the centre of the bin. For the red curve, the exact k-value 
for each mode is used when computing the sum in each bin. 
For small bins, or fields where the values vary slowly as a 
function of k, these two should agree and indeed they do 
for smaller scales where our (logarithmic) bins are smaller. 
There is a difference between the methods for the largest 
scales in our simulations, this will be discussed below for 
each test. 

For the density field, the two methods for calculating 
the ratio do give different answers. However, for both meth¬ 
ods, the deviation is within 2% for every bin except the first. 
Thus, this consistency check for the density field is well sat¬ 
isfied for all scales k ^ 0.2/iMpc -1 . 

The consistency check for the velocity field is less well 
satisfied: there is a sharp divergence in the power spectra 
on the smallest scales, such that the check is not satisfied 
within 10% at k ~ 8/i _1 Mpc. This shows the effect of the 
DTFE window function on the extracted fields. We will not 
consider the extracted vector potential for k larger than 
k ~ 8/i -1 Mpc when presenting our results. Furthermore, the 
two methods show very different behaviour: the method us¬ 
ing the average fc-value for each bin causes the consistency 
test to fail on large scales. However, with the more exact 
method, the consistency check is very well satisfied on all of 
the largest scales. This suggests that the dominant contri¬ 
bution to the bins on the largest scales comes from the low-A: 
end of each bin, hence the overestimation of k 2 P$(k) when 
the average fc-value for each bin is used. The strong effect 
here is partly caused by the relatively steep slope of the ve¬ 
locity power spectrum. We note that this effect would also 
come into play when calculating the dimensionless veloc¬ 
ity power spectrum for binned data. Nonetheless, the good 
agreement of the consistency check when using the second 
method is strong evidence that the derivatives of the velocty 
field are being calculated correctly. 

A further check that we can perform is to extract the 
complete momentum field, p = (1 + S)v, and decompose it 
into its vector and scalar parts directly rather than dealing 
with derivatives. The power spectrum of the vector potential 
can then be calculated from the vector part of the momen¬ 
tum field, p", using 

p SN (k)=( i6 ^rJ pMk). (i4) 

In figure [T] we show the ratio of the vector power spectrum 
calculated using the two methods, with the different lines 
corresponding to different individual realisations. The vector 
potentials calculated from the two methods are broadly con¬ 
sistent, within 20% for most of the range under considera¬ 
tion, and agreeing to within a factor of 2 for k ^ 0.2/iMpc^ 1 . 
We are unsure what the causes of the difference between the 
two methods are. In particular, we checked for whether there 
is an effect coming from the use of k averaged over the bin, as 
in the velocity field consistency check, however this effect is 
negligible for the gravitomagnetic potentiaQ The difference 

1 As an aside, we note that we also calculated the momentum 
field by extracting the velocity field and density field separately 
at each grid point, before multiplying them together. The power 
spectrum calculated from this field agrees well with that calcu- 
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between the methods is larger than the variation amongst 
realisations for either method. 

We can also extract the momentum held directly us¬ 
ing a standard cloud-in-cells (CiC) approach ( |Hockney fe| 
Eastwood 1981), and compare this to the momentum held 


extracted using the DTFE code. The ratio between these 
helds is shown in figure [2] There is good agreement between 
the two methods of computing the momentum power spec¬ 
trum on larger scales, but with a divergence between the two 
methods on smaller scales. It is unclear which method would 
be expected to be more accurate on these smaller scales: 
the DTFE method suffers from having a window function 
that cannot be deconvolved, however the CiC method will 
have cells with a zero momentum held, due to the lack of 
nearby particles, for a sufficiently large grid. In fact, the CiC 
method does not converge as the grid size is increased. We 
used a 512 3 grid for the CiC code, although we checked that 
changing this to 256 or 1024 does not signihcantly affect the 
results. Unlike the DTFE method, derivatives cannot be di¬ 
rectly extracted with the CiC method, so the consistency 
checks performed earlier for the DTFE method cannot be 
applied to the CiC method. This also means that the first 
method of extracting the vector potential, using the curl of 
the momentum field, cannot be carried out with the CiC 
method. 

We present the vector power spectrum from both the 
momentum field and the curl method in the results section. 
We note that the level of agreement between figures]]] and [2] 
suggests that our vector potential power spectrum is robust 
and correct to within a factor of 2. It is unclear to us which 
method should be trusted more; whilst the momentum held 
method is simpler, the derivative method allows us to ex¬ 
amine the different components, notably the vorticity, and 
check that it behaves as expected. The differences between 
the two methods do not affect the observability of the vector 
potential, see Bruni et al. (20141 and Thomas et al. (20141. 


2.2.5 Comparison to previous findings 


There are several works in the literature to which we can 
compare our Endings on the velocity held and its compo¬ 
nents. As mentioned above, the vorticity and velocity di¬ 
vergence power spectra were extracted from N-body simula¬ 


tions in Pueblas & Scoccimarro (20091 using an alternative 


implementation of the Delauney tesselation. They found a 
strong dependence on resolution of the extracted vorticity 
power spectrum and an approximate scaling of the vorticity 
power spectrum with the seventh power of the linear growth 
factor. 

The vorticity and velocity divergence power spectra in 


Pueblas & Scoccimarro (2009) are consistent with the spec¬ 


tra extracted for this paper and we found the same res¬ 
olution dependence of the vorticity power spectrum (see 


above). However, as detailed in Appendix B2 we do not 


End the same scaling of the vorticity spectrum with the sev¬ 
enth power of the linear growth factor (D+): Although this 


Gravitational potential field: momentum field vs derivatives 



Figure 1. The ratio of the vector potential power spectra com¬ 
puted using the vector part of the momentum field and the curl 
of the momentum field. The blue, magenta and red curves show 
the ratio for the three realisations of the HR simulations, and the 
black (dashed) curve shows the average over these three. There 
is reasonable agreement between the two power spectra for the 
smaller scales, however the two methods diverge for the largest 
scales and there is a difference of a factor of 5 at the largest scales. 
For most of the range of k under consideration (k ^ 0.2fiMpc —1 ), 
the two vector power spectra agree to within a factor of 2. 


Vector part of momentum power spectrum: CiC/DTFE 



Figure 2. The ratio of the vector potential power spectra com¬ 
puted using the vector part of the momentum field calculated us¬ 
ing the Cloud-in-Cells method and the DTFE method. The blue, 
magenta and red curves show the ratio for the three realisations 
of the HR simulations, and the black curve shows the average over 
these three. The two methods agree very well on larger scales, but 
diverge for the smallest scales. 


lated by extracting the momentum field as a single field. The 
same agreement is not obtained when extracting the field S 2 and 
comparing to squaring the density field, when using either the 
DTFE code or a CiC method. 
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scaling seems to hold at low redshift, it no longer holds at 
redshift one and beyond. At these earlier times, the power 
spectrum is smaller than expected from the growth factor 
to the seven scaling, so the vorticity power spectrum must 
have grown by less at redshift two than expected. 

Two recent publications (Zheng et al.|2013 |Koda et al. 


|2013| ) have examined the velocity held from the point of 
view of redshift space distortions. In these works, a differ¬ 
ent method of extracting velocity fields is used, the nearest 
particle method. In this method, the velocity at a grid point 
is given by the velocity of the nearest particle to that grid 
point. See those works for comments on the differences be¬ 
tween the nearest particle and Delauney tesselation methods 
of extracting the velocity power spectra. Here, we note that 
there appear to be pros and cons to both methods, with 
no clear “better” method. It would be interesting to exam¬ 
ine how close the agreement between the vector potentials 
extracted by the DTFE and nearest particle methods is. 

Nonetheless, there are some general observations that 
can be compared between these works. Notably, the magni¬ 
tude of the velocity and vorticity spectra is found to be sim¬ 
ilar, considering the differences in cosmological parameters. 
Also, the onset of non-linearity is found to occur at lower 
k for the velocity divergence than for the density. In addi¬ 
tion, Zh eng et al.j (J2013) finds a strong dependence of the 
curl component of the velocity field on the resolution, simi¬ 


larly to both this paper and Pueblas & Scoccimarro (20091. 


They also find a time dependence of this component that is 
approximately D‘ + up to z = 2, although this relationship 
breaks down by up to a factor of two for certain redshifts 
and scales. As mentioned above, whilst our simulations also 
find this time dependence of the vorticity at low redshift, we 
find that the relationship breaks down for z > 1. There is no 


examination of multiple realisations in Zheng et al. (2013) 
and, similarly to the comments in Appendix |B2| regarding 


Pueblas & Scoccimarro (2009), the difference between our 


realisations is sufficient to explain the difference between our 


results and those of Zheng et al. (20131. 


The broad agreement between different methods, in¬ 
cluding agreement regarding resolution dependence and con¬ 
vergence, is promising. Details of the vorticity field and its 
evolution require further study, but the vorticity is a sub¬ 
dominant contribution to the vector potential. As the simu¬ 
lations and snapshots used in the papers mentioned in this 
section are different to ours, it is not possible to compare 
the methods and extracted fields any more precisely. We 
note that the three works mentioned here do not have mul¬ 
tiple realisations of their high resolution simulations, so we 
are unable to determine if the variation in vorticity between 


realisations found by us is reproduced (see Appendix B4 1 . 


As this manuscript was being prepared, |Hahn et al.| 
|2014| appeared on the arxiv. This paper investigates the 
properties of velocity divergence and vorticity and confirms 


many of the findings of Pueblas & Scoccimarro (2009). In 


particular, they agree with our results regarding the conver¬ 
gence of the DTFE code for sufficient mass resolution and 
our finding of a resolution dependence of the velocity di¬ 
vergence, which did not appear in |Pueblas fe Scoccimarro| 
12009). They use a different method to compute the vortic¬ 


ity and velocity divergence power spectra, which agrees with 
the DTFE code for sufficient resolution. However, as with 
the previous papers, there seems to be no examination of 


multiple realisations with the same resolution, in order to 
compare our findings. In addition, there is no examination 
of the time dependence and thus no confirmation or rejec¬ 
tion of the D\ scaling of the vorticity spectrum at higher 
redshifts. 


3 RESULTS 

In this section we present the power spectrum of the post- 
Friedmann vector potential as calculated from N-body sim¬ 
ulations. We show the power spectrum at z = 0 and the 
different components of the source, as well as the evolution 
of the power spectrum between z = 2 and z = 0. In addi¬ 
tion, we show the ratio between the vector and scalar power 
spectra, and examine the time evolution of this quantity 
as well. The power spectra plotted for the scalar and vector 
gravitational potentials are the dimensionless power spectra. 
The closest analytic result to our calculation is the second 
order perturbative vector potential calculated in |Hui-Ching| 
Lu et al. (2009). We will compare our results to theirs at 


redshift z = 0, as well as comparing the time evolution. 


3.1 Results at redshift zero 

In figures [3] and [4j we show the power spectrum of the post- 
Friedmann vector potential as well as the standard Newto¬ 
nian scalar potential, at z = 0, for the curl and momen¬ 
tum field methods of extraction respectively. As expected, 
both methods show that the scalar potential is small over 
all scales and the vector potential is subdominant. There is 
a quantitative difference between the two methods on the 
largest scales, but this difference is not sufficient to alter the 
expected qualitative behaviour. Notably, the effect of the 
vector potential on weak-lensing power spectra, as examined 


Thomas et al. (2014), will remain negligible, regardless of 


which method is used to calculate the vector potential. We 
have been unable to determine the reason for this discrep¬ 
ancy and it is unclear to us which method should, a priori, 
be expected to be more accurate. 

In figures [5] and |6j we show the ratio between the power 
spectra of the vector and scalar gravitational potentials at 
redshift zero, for the two methods of extracting the vec¬ 
tor potential. We plot the ratios for all three individual re¬ 
alisations of the HR simulations. For the curl method, as 
shown in Bruni et al. pollt , this ratio is approximately 
2.5 x 10 _s . This ratio does not vary significantly over the 
range of scales considered, although there is a slight increase 
towards smaller scales. However, for the momentum field 
method, the ratio is not approximately constant due to the 
decreased power on large scales. We will compare this be¬ 
haviour to the analytic second-order perturbative behaviour 
shortly, here we just note that the curl method produces 
qualitative behaviour that is closer to the analytic predic¬ 
tion. 

In figure [7J we show the power spectra of the three 
sources of the vector potential using the curl method, see 
equation |8|. The power spectra plotted here are given by 
P(k)/ (/ 2 7d 2 (2-7 t) 3 ), where T-L is the conformal time Hubble 
constant and / = din D /din a is the logarithmic derivative 
of the linear growth factor D. These units are chosen such 
that the power spectrum of the velocity divergence agrees 
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Gravitational Potentials 



Figure 3. The scalar (dashed red line) and vector (solid blue 
line) gravitational potential power spectra at redshift zero, with 
the vector potential calculated using the curl method. The lin¬ 
ear theory scalar potential is shown for comparison (dotted black 
line). 


Ratio of Gravitational Potentials 



Figure 5. The ratio at redshift zero between the vector potential, 
calculated using the curl method, and the scalar potential. The 
three curves show the ratio for the three realisations of the HR 
simulations. 


Ratio of Gravitational Potentials 



Gravitational Potentials 



k/h Mpc -1 


Figure 4. The scalar (dashed red line) and vector (solid blue 
line) gravitational potential power spectra at redshift zero, with 
the vector potential calculated using the momentum field method. 
The linear theory scalar potential is shown for comparison (dotted 
black line). 


with the density power spectrum on linear scales and have 
the same units as the matter power spectrum, following 


Pueblas & Scoccimarro (2009). The vorticity, although often 


ignored in perturbation theory, is the only one of these three 
quantities that is linear in perturbations. This figure shows 
that it is negligible compared to the other two components, 
so the vector potential is being predominantly generated by 
non-linear effects. 

Since this vector potential is the first correction to New¬ 
tonian theory, this calculation is the first quantitative check 


Figure 6. The ratio at redshift zero between the vector potential, 
calculated using the momentum field method, and the scalar po¬ 
tential. The three curves show the ratio for the three realisations 
of the HR simulations. 


of the relationship between Newtonian simulations and GR 
on fully non-linear scales. The small magnitude of the vec¬ 
tor potential suggests that running Newtonian simulations 
is sufficiently accurate for cosmological purposes, whereas a 
larger calculated value for the vector potential would sug¬ 
gest that the approximations taken in deriving the fully non¬ 
linear Newtonian equations do not hold sufficiently well. As 
far as relating Newtonian and relativistic cosmologies goes, 
in the language of Green & Wald (2012), the smallness of 
this vector potential allows the use of the abridged dictio¬ 
nary in |Chisari &; Zald arriaga (2011), rather than the dic¬ 
tionary proposed in Green & Wald (2012). We note that the 
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Power spectra for sources of vector potential Redshift Evolution of vector potential 




Figure 7. The power spectra of the three source terms for the 
vector potential in equation [8] , the vorticity (dashed black line), 
V<5 X v (dot-dashed blue line) and 5V X v (solid red line). The 
power spectra plotted here are given by P(k)/ (2n) 3 ), such 

that the power spectrum of the velocity divergence agrees with 
the density power spectrum on linear scales and ensuring that all 
of the power spectra have the same units, following Pu eblas fe| 
|Scocci marro ( 2009]) . The linear matter power spectrum is shown 
as a dotted magenta line for comparison. 


analysis here is for a ACDM cosmology, further work is re¬ 
quired to determine the validity of Newtonian simulations 
in general dark energy cosmologies. 


3.2 Time evolution 

In this section we will examine the time evolution of the 
vector potential, and its ratio to the scalar potential, for the 
redshifts listed in table[2] The vector potential is this section 
has been computed using the curl method. In figure |8j we 
plot the ratio of the vector potential to the scalar potential 
as a function of redshift. The different curves in this plot 
show the evolution for different wavenumbers. We can see 
that individual k-modes do not exhibit significant growth 
over time, although the more non-linear scales do exhibit 
slightly more variation in time. Similarly to the scalar grav¬ 
itational potential, the vector potential at a fixed scale is 
not monotonic over time on non-linear scales. 

In figure [9] we plot the ratio of the vector potential to 
the scalar potential as a function of redshift. The different 
curves in this plot show the same wavenumbers as in figure 
[5] The ratio stays fairly constant over time, varying by less 
than a factor of two for a given scale. Across the entire range 
of times and scales under consideration, the ratio varies by 
less than a factor of 4. The ratio between the gravitational 
potentials is also not monotonic over the redshift range un¬ 
der consideration for a given scale. 


3.2.1 Comparison to perturbative calculation 

In |Hui-C hing Lu et ah] ([2009]), an analytic calculation of the 
vector potential was performed using perturbation theory. 


Figure 8. The evolution of the vector potential for six different 
wavenumbers. From top to bottom, these are k = 0.23hMpc —1 
(brown), k = 0.55fiMpc —1 (black), k = 0.79/iMpc —1 (cyan), 
k = 1.01/iMpc -1 (blue), k = 2.51/iMpc -1 (magenta) and k = 
5.03/iMpc -1 (red). 


Redshift Evolution of Gravitational Potential ratio 



Figure 9. The ratio of the vector potential to the scalar po¬ 
tential plotted for six different wavenumbers. From bottom to 
top (at redshift=l), these are k = 0.23/rMpc" 1 (brown), k = 
0.55/iMpc -1 (black), k = 0.79hMpc -1 (cyan), k = l.OlhMpc -1 
(blue), k = 2.51/iMpc —1 (magenta) and k = 5.03hMpc —1 (red). 


As a perturbative analysis, it is unclear how large a value 
of k this calculation should be extended to. Here we will 
assume it is valid on all of the scales of overlap between this 
method and ours. 

For the curl method of computing the vector power 
spectrum, there is similar qualitative behaviour between the 
two methods, with the ratio of the power spectra of the vec¬ 
tor and scalar potentials being fairly constant and of or¬ 
der 10“ 5 in both methods. The difference between the two 
methods being that the ratio in Hui-Ching Lu et al. (20091 
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is between the vector and the linear theory scalar potential, 
whereas the our ratio is between the vector and the fully 
non-linear scalar potential. This means that despite this sim¬ 
ilar qualitative behaviour, the power spectrum of the vector 


potential in Hui-Ching Lu et al. (2009) underestimates the 


fully non-linear value on these scales by up to two orders 
of magnitude, the same factor by which the linear theory 
scalar potential power spectrum underestimates the power 
spectrum of the fully non-linear scalar potential. 

The momentum field method of calculating the vector 
power spectrum results in less-similar qualitative behaviour. 
It is unclear how well the gravitomagnetic potential would 
be expected to match the perturbative prediction on these 
scales as the velocity field differs from the linear theory at 
larger scales than the density. 

The power spectrum of the perturbative vector poten¬ 


tial is given in Hui-Ching Lu et al. (2009) as 


Vs(k) = 


( 2 A tc \ 4 (3g W+U}] 


V 5 S° 


9raW 2 


fc 2 n(u 2 ) 


(15) 


where V s is the dimensionless power spectrum of the vec¬ 
tor potential, A 7 z is the primordial power of the curvature 
perturbation, g is the growth factor for the scalar potential, 
goo is a normalisation parameter chosen so g( 0 ) = 1 , n is a 
function of the transfer function, is the time dependent 
matter density and 'H is the conformal Hubble constant. The 
second term in parentheses contains all of the time depen¬ 
dence of the vector potential power spectrum and essentially 
acts as the growth factor for the vector potential. We have 
compared this perturbative prediction for the growth factor 
of the vector potential to the growth measured in the simu¬ 
lations (see figure 5 in file CC). This shows that the analytic 
prediction is not the main source of the time evolution of the 
vector potential. 


4 CONCLUSION AND DISCUSSION 


In this paper we have presented the post-Friedmann 
frame-dragging vector potential calculated on non-linear 
scales from N-body simulations. We have presented this 
vector potential at redshift zero, as well as examining its 
evolution with redshift. We have also presented the tests 
we have performed in order to establish the robustness of 
our result, including tests of simulation parameters and 
different methods of extracting the source of the vector 
potential. 

We have shown that our density, velocity divergence and 
vorticity spectra are consistent with the literature and show 
similar behaviour regarding convergence tests, particularly 
mass resolution. We do not see the vorticity scaling with 
the seventh power of the linear growth factor D+ ( |Pueblas] 
fe Scoccimarro|200 9) beyond 2 = 1, however the differences 
between our results and others’ are within the variance 
between realisations. We have noted a larger variation of 
the vorticity than the density and velocity divergence fields 
between different realisations, a result that does not seem 
to have been studied in the literature. 

We have shown that there is no evidence for a systematic 
dependence of the vector potential spectrum on box size for 
boxes smaller than 200/i - 1 Mpc, or on mass resolution with 
mass resolution better than 6 x 1 O 8 M 0 . There is also no 


evidence that the vector potential is sensitive to the soften¬ 
ing length, binning, number of samples (an internal DTFE 
parameter) or the grid size used in the analysis. There is a 
reasonable agreement between the different methods (curl 
and momentum field) of extracting the vector potential, 
although there is an unresolved discrepancy between the 
two methods on the largest scales. We do however note the 
importance of the variation of the vector potential between 
realisations, this issue is discussed more fully in Appendix 


Figures [3] and [ 8 ] comprise the main physical results of 
this paper, showing the magnitude of the vector potential 
power spectrum at redshift zero and its evolution with time 
respectively. The magnitude of the vector potential power 
spectrum can also be expressed in terms of its ratio to 
the power spectrum of the scalar potential, as shown in 
figures [5] and [9] We have shown that the power spectrum 
of the vector potential is around 10 5 times smaller than 
the power spectrum of the scalar potential, over a range 
of scales and redshifts. These values were used in IBrunil 


et al. (20141 and Thomas et al. (20141 when examining 


the observability of the vector potential, showing that it 
is neglgible for currently planned weak-lensing surveys. 
The small magnitude of the vector potential found here is 
the first quantitative check of the validity of Newtonian 
simulations compared to GR on fully non-linear scales and 
supports the use of Newtonian simulations for computing 
cosmological observables. In terms of interpreting the 
simulations, the small value of this vector potential seems 
to justify the use of the abridged dictionary in Chisari fe| 


Zaldarriaga (20111, rather than the dictionary proposed 


Green & Wald (2012 1 , for relating GR and Newtonain 


cosmologies. 

The work carried out so far considers a ACDM cosmology, 
so this conclusion may no longer be true for a dark energy or 
modified gravity cosmology. The post-Friedmann approach 
would need to be expanded to include modified Einstein 
equations and/or a fluid with pressure in order to examine 
alternative cosmologies and determine whether the use 
of Newtonian-type N-body simulations is still valid in 
those cosmologies. The post-Friedmann expansion has been 
applied to f(R) gravity and the vector potential calculated 
from f(R) simulations in Thomas et al. (20151. The vector 
potential in f(R) was found to be larger than in General 
Relativity. We hope that this, and further extensions to 
the work in this paper, will allow us to understand how 
generic the findings in this paper are, and thus justify 
one of the most widely used tools in cosmology, N-body 
simulations. Whilst this manuscript was being prepared for 
submission, Adamek et al. (20141 appeared on the ArXiv. 
Their preliminary results seem to agree with the results of 
this work. It will be interesting to perform a more in-depth 
comparison once the details of their work are available. 
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APPENDIX A: VECTOR POWER SPECTRA 

We will be dealing with vector quantities, for which there 
are different ways to define the power spectrum. Our power 
spectrum for a generic vector v is defined as 

(v(k) • v (k 1 )) = (2n) 3 S 3 (k — k')Ps(k) (Al) 

Note that for a divergenceless vector, such as B N , 
k 2 Pg N (k) = P V xfl JV (^)- With our Fourier transform con¬ 
vention, the dimensionless power spectrum for a field X is 
given by Ax = k 3 Px(k)/2n 2 . All plots of the power spec¬ 
trum of the vector potential show the dimensionless power 
spectrum. 

Using equation 0 

(V x B N (k) ■ V x B N *(k')) = ( ^ 67V ^ ba " S j 

([(V5)~x F+(l +7fvx«] • [(V^J+fl + ijvxt?]*) (A2) 

(vxi f^(k) ■vxB N *(k , ))=( i67r ^ 2 pi,a2 y 

f ([VS>Cv\ • [vTxv\) +([v7xV| • \8V~xv\) -t-<[V^r>TV| • [VxFv\) N 
+{[8sTxv I ■ (V^xv\) +([SV>Cv\ • [txTxvf) -t-<[AV^TV] • [VxAv\) 

yf([Vx^] ■ [v7xV|*> + ([VxFv\ • [sVx'v ]*) + ([V~xTi>] • [VxFv\ ) 

Noting that A ■ B* = {A* ■ B)*, 

([V~xii] • [vTxv\) + ([vLTw] ■ [Vxi/]}* = 

2re (((V x v) • (V<S x v)}) = (2?r) 3 5‘ ! (k - k 7 )P (V 5xv)(vxV)(kXA3) 

And therefore the dimensionless power spectrum for the 
vector potential is given by 

A Mk) = ( 16 tr 2 ) 2 ^^(^ (A4) 

where 

Pdv(k) — Pv xv(k^ T p5vx^(^) T Pcvs)xfi(k) 

+P(VSx »)(VXJ)W + P(VSXv)(SVXv)(k) + P(SVX v)(VXv)(k) (A5) 


APPENDIX B: ADDITIONAL ROBUSTNESS 
INFORMATION 

In this Appendix we show the figures referred to in the main 
text as well as discussing additional robustness and conver¬ 
gence tests that were carried out in order to establish our 
result. 


Resolution Dependence: Vorticity 



Figure Bl. The vorticity power spectra extracted from simu¬ 
lations with varying box size and mass resolution. Lines with 
the same colour share the same mass resolution (in units of 
10 8 M 0 : 3.97 (red), 6.26 (magenta), 6.31 (yellow), 7.76 (green), 
16.3 (cyan), and 130 (blue). The black curve is the linear matter 
power spectrum for comparison. 


Resolution Dependence: Vector Potential 



Figure B2. The vector potential power spectra extracted from 
simulations with varying box size and mass resolution. Lines 
with the same colour share the same mass resolution (in units 
of 10 8 Mq: 3.97 (red), 6.26 (magenta), 6.31 (yellow), 7.76 (green), 
16.3 (cyan), and 130 (blue). The spectra have been multiplied by 
k 2 in order to better show the variation. 

Bl DTFE parameters 

There are several internal DTFE parameters that are used 
when computing these fields on a regular grid. We investi¬ 
gate the effects of two of these parameters here, the grid size 
and the number of samples that are made in each grid cell, 

N samples • 

We examined the effect of varying the grid size on the 
extracted density, velocity divergence and vorticity power 
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Box Size Dependence: Vector Potential 



k/h Mpc -1 

Figure B3. The vector potential power spectra extracted from 
simulations with varying box size and mass resolution. Lines 
with the same colour share the same box size: 140/i -1 Mpc 
(red), 160/i _1 Mpc (magenta), 200/i _1 Mpc (yellow), 240h -1 Mpc 
(green), 320h~ x Mpc (blue). The spectra have been multiplied by 
k 2 in order to better show the variation. 


Ratio of 100 samples to 1000 samples 



Figure B4. The ratios of the power spectra computed with 
^samples = 100 and IVsampies = 1000. The ratios shown are for 
the density (red), velocity divergence (blue), vorticity (cyan) and 
vector potential (black). 


spectra. In all cases the agreement is very good, except 
on the smallest scales. A discrepancy on this scales is ex¬ 
pected due to the change in the resolution of the grid and 
the effects of the DTFE window function. However, even 
on the smallest scales, the discrepancy is small. This can 
be seen in figures 1, 2 and 3 in file BG, where we show 
the extracted spectra from one of the 160h _1 Mpc simu¬ 
lations with 1024 3 particles at redshift zero. The different 
lines show the different grid sizes used: 1024 (blue), 950 
(cyan), 850 (green), 750 (magenta) and 640 (red). The black 
line shows the linear density power spectrum for compari¬ 
son. For our results plots, we have used the suggested value 
grid = Alpart = 1024. 

Our analysis has all been carried out with Al samp i es = 
100 points per grid cell, partly due to computing constraints; 
increasing the number of samples increases the run time and 
memory required when analysing a snapshot. However, in 
figure |B4| we show the effect of increasing Al samp i es to 1000 
points per grid cell for one of the 160h _1 Mpc simulations 
with 1024 3 particles. The velocity divergence and vorticity 
spectra agree very well between the two different numbers 
of samples. The density power spectrum shows a deviation 
that increases towards smaller scales, however is within 5% 
for the range of scales under consideration here. The power 
spectrum of the vector potential shows more deviation, with 
decreasing deviation for smaller scales. However, the change 
in the vector potential is within 10% for every bin after the 
first and is within 5% for all scales k > 0.3/i _1 Mpc. 


B2 Linear evolution 

A further check that can be performed is to examine how the 
time variation of our extracted density, velocity divergence 
and vorticity power spectra compares to the respective linear 
predictions. For the density and velocity divergence fields, 


the power spectra evolve as (D+(z)/D+(z = 0)) 2 on the 
largest scales and earliest times, as per the linear theory 
prediction. This prediction becomes increasingly inaccurate 
for more scales at lower redshifts due to non-linear effects. 

The time evolution that we found for the vorticity is 
shown in figure |B5| In this figure, the power spectrum at 
each redshift has been divided by the seventh power of the 
linear growth factor for that redshift, (D + (z)/D + (z = 0)) 7 , 
as suggested by Pueblas fe Scoccimarro|([2009||. In|Pueblas fe 


Scoccimarro (2009), the authors include an approximate an¬ 


alytic derivation of the time evolution of the vorticity power 
spectrum, finding it to behave as f4(z)D+(z)), where fv(z) 
is the fraction of the volume that undergoes orbit crossing. 
Fitting to their simulations, they found (D + (z)/D + (z = 


0)) 7 


to be the best fit value. The scaling of our vortic¬ 


ity spectrum appears similar to that found in |Pueblas fe| 
|Scoccimarro| ( |2009[ ) . However, in our simulations this scal¬ 
ing appears to break down for higher redshift, 2^1. We see 
a smaller vorticity spectrum at these times than expected 
from the (D+(z)/D+(z = 0))' scaling. Figure 


B5 


shows this 

discrepancy along with the error amongst our simulations. 
These errors do not appear sufficiently large to explain the 
discrepancy. However, it is worth noting that the variation 
amongst our realisations (see Appendix B41 is large enough 
to explain the difference in the time evolution of the vortic¬ 
ity between our simulations and the single high resolution 
simulation in Pueblas & Scoccimarro ( 2009). The variation 
between realisations was not considered in IPueblas fc Scoc-I 
[cimarro ] I?!? 09 ), however it seems likely that the function 
fv(z) varies between realisations. The range of the scaling 
of the vorticity with the linear growth factor has an upper 
value of 7.3 in P ueblas fe Scoccimarro| (2009). Using this 
value reduces, but does not remove the discrepancy. 
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Time Evolution: Vorticity 



Figure B5. The non-linear vorticity power spectrum at selected 
redshifts, z = 2 to z = 0, each divided by the respective linear 
theory density growth factor to the seventh power, see |Pueblas| 
|fe Scoccimarro| | |2009| l. See table [ 2 ] for details and an explanation 
of the colours. As expected, the linear theory prediction works 
well on the largest scales and is generally worse for smaller scales 
and later times, however the scaling as the seventh power of the 
density growth factor seems to break down at earlier times. The 
error bars on this plot show the standard error on the mean for 
each set of realisations. 


B3 Comparison with the POWMES density 
power spectrum 

The density and density gradient power spectra (the lat¬ 
ter divided by fc 2 , see the consistency check above) that 
we have extracted can be compared to the density held 


extracted by POWMES (Colombi et al. 2009), a state of 


the art conventional density power spectrum estimator. For 
the HR simulations, the power spectra agree within 10% for 
0.2/iMpc -1 ^ k ^ 7.0hMpc _1 , see figure 3 in hie CC, and 
within 5% for the majority of this range. A similar result is 
seen for the ratio of the DTFE gradient of the density spec¬ 
trum (divided by fc 2 ) to the POWMES density spectrum, 
see hgure 4 in hie CC. 

The agreement on the largest scales, in the hrst 4-5 bins, 
is affected by the choice of binning. If the number of bins 
used for the DTFE extraction is doubled, then the DTFE 
and POWMES extractions agree much more closely as the 
bins are then of a more similar size and location. As noted 
in Appendix |B7[ if we increase the number of bins then 
the number of k modes contributing to the hrst few bins is 
much smaller, so we will continue to use 7V gr ia/4 bins in our 
analysis. The agreement between the POWMES and DTFE 
methods is sufficient to support the robustness of our density 
and density gradient spectra. 


B4 Realisations 

In this section we show how the extracted spectra vary 
amongst realisations. We will illustrate this with the 
160h _1 Mpc 640 3 particle simulations for which there are 8 


realisations. In all cases we consider the variation at redshift 
zero. 

We examined the variation amongst realisations for the 
density held, using both the DTFE code and POWMES, and 
also the velocity divergence. These all showed the expected 
variation, namely that cosmic variance causes a difference 
between the realisations on the largest scales in each box, 
but this difference is much reduced on smaller scales. The 
variance between realisations for the density held was very 
similar for the two methods of extracting the density held. 

In hgure |B6| we show the variation of the vorticity held 
amongst realisations. This plot shows that the variation 
amongst realisations is greater for the vorticity than for the 
density. On smaller scales, the variation amongst realisations 
of the vorticity is less than on large scales, but still greater 
than for the density Held. We are not aware of this being 


previously noted in the literature, and the works (Pueblas 


& Scoccimarro 2009 Zheng et al. 2013 Hahn et al. 2014 1 


that we compare our vorticity spectrum to in the main text 
do not have multiple realisations in order to have seen this 
effect. 

In hgure |B7| we show the variation amongst realisa¬ 
tions of the vector potential. On large scales, the variation 
between realisations is very similar to that between the vor¬ 
ticity spectra. However, the variation does not appear to re¬ 


duce on smaller scales. According to perturbative results (Lu 


et al. 2008 Hui-Ching Lu et al. 2009), the vector is generated 


most efficiently by coupling between two different fc modes, 
particularly if one of them is entering the horizon. Given the 
similar qualitative behaviour of the fully non-linear vector 
potential, it is reasonable to assume that this is also gener¬ 
ated by coupling between large scale modes and small scale 
modes. Thus, the large scale variance between realisations 
will be affecting the vector power spectrum on smaller scales, 
resulting in the variance between realisations not decreasing 
on small scales. 

In hgure [B8l we show how the value of the vector poten¬ 
tial from the individual realisations of the HR simulations 
compares to the average over realisations of simulations with 
different parameters. Note that the variation between the 
HR realisations is greater than the variation between the 
average over realisations for different simulation parameters. 

As mentioned above, the increased variance between re¬ 
alisations may be an unavoidable feature of the vector po¬ 
tential. As such, this represents the dominant source of error 
in calculating the vector potential, as long as the simulation 
parameters are sufficiently good. If an observational test of 
the vector potential was found, then many more realisations 
than the number carried out for this paper would be re¬ 
quired, in order to more carefully investigate this effect and 
determine more precisely what the observational prediction 
would be for a ACDM cosmology. 


B5 Softening length 

In this paper we have chosen our softening lengths follow¬ 
ing 


Pueblas & Scoccimarro (2009) in order to compare to 


their results. In figure [B9l we show how a 160/i -1 Mpc sim¬ 
ulation with 640 3 particles and the same initial conditions 
varies if the softening length changes from 6.5kpc to 5kpc. 
This is a 20% change in the softening length. The variation 
between the density, velocity divergence and vorticity spec- 
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Vorticity Realisations Vector Potential Convergence 




Figure B6. The vorticity power spectra as extracted from the 
8 160h —1 Mpc simulations with 640 3 particles. The dashed black 
line denotes the average of the 8 simulations. 


Vector Potential Realisations 



k/ h Mpc 1 

Figure BT. The vector potential power spectra as extracted from 
the 8 160h —1 Mpc simulations with 640 3 particles. The dashed 
black line denotes the average of the 8 simulations. The spectra 
have been multiplied by k 2 in order to better show the variation. 

tra is very small for this change. The power spectrum of the 
vector potential varies more, but is within 5% of the value 
for nearly the entire range under consideration. Since this 
5% variation is significantly smaller than the 20% variation 
in the softening length, we do not think the choice of soft¬ 
ening length significantly impacts our results for a sensible 
choice of softening length. 

B6 Smaller Boxes 

Here we examine some additional plots that demonstrate 
some of the comments made in the main text. We ran a 
set of 8 simulations with 512 3 particles in an 80h -1 Mpc 


Figure B8. The vector potential power spectra from different 
simulations, divided by the average vector potential from the 
three 160/i —1 Mpc simulations with 1024 3 particles. The three red 
curves show the vector potential from the three realisations of 
the 160/i — x Mpc simulations with 1024 3 particles. The cyan and 
magenta curves show the vector potential from the average of the 
160h~ 1 Mpc simulations with 640 3 and 880 3 particles respectively. 
The yellow curve shows the average of the 140h —1 Mpc simula¬ 
tions with 768 3 particles and the green curves shows the average 
of the 200h —x Mpc simulations with 1024 3 particles. Note that 
the variation between the high resolution simulations is greater 
than the variation between the average values from simulations 
with different parameters. 


Varying Softening 



k/ h Mpc -1 


Figure B9. The ratio between the power spectra extracted at 
redshift zero for the same initial conditions run with two differ¬ 
ent softening lengths. The different spectra plotted are density 
(red), velocity divergence (blue), vorticity (cyan) and the vector 
potential (black). 
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box and extracted the power spectra in the same way as 
from our other simulations. In figures 3, 5 and 6 in file RB, 
we show the density, vorticity and vector potential power 
spectra respectively, colour coded to match box size as in 
figure [B3] In addition, the power spectra extracted from the 
smaller boxes is shown as a dashed black line. It is clear that 
the spectra extracted from the smaller 80h _1 Mpc boxes are 
systematically smaller, irrespective of any other dependence 
on box size and resolution. The effect of using (too) small 
boxes when running N-body simulations has been examined 


in the past, see e.g. Gelb & Bertschinger (1994); Bagla & 


Ray (2005 I; Bagla & Prasad (2006). A suggestion in Bagla & 


Prasad (20061 is that it is important that the ratio between 


the box size and the scale of non-linearity is sufficiently large. 
As a result, for our simulations, the smallest boxes we have 
run that we consider trustworthy are the 140/i^ 1 Mpc simu¬ 
lations. It remains to be seen wheth er smaller boxes, suc h as 
the 100h _1 Mpc simulations used in Zheng et al. |2013l are 
a robust source of spectra such as the density and vorticity. 


B7 Number of bins 

We considered the effect on our extracted power spectra 
of varying the number of bins. As expected, increasing the 
number of bins increases the noise of the power spectra and 
there is no systematic deviation. Our results from varying 
the number of bins on the density, velocity divergence and 
vorticity power spectra are shown in figures 4-6 in hie GB. In 
each of these plots, the 256 bins used for the analysis in this 
paper is shown by the black line, the blue lines denote the 
use of 512 bins and the red lines are for 1024 bins. We have 
used 256 bins for our analysis to ensure that the low k bins 
contain a sufficient number of fc-modes. For the 256 bins, 
the first two k bins contain 58 and 218 fc-modes respectively, 
whereas these numbers are 12 and 41 for the corresponding 
bins when 1024 bins are used. Note that, as mentioned in the 
POWMES section, the variation between the 256 bin and 
512 bin power spectra is similar to the variation between 
the POWMES method and the DTFE method using 256 
bins. This is due to the number and location of bins in the 
POWMES method being very similar to the DTFE method 
with 512 bins. 

In addition, figure 7 in hie GB shows the variation of 
the vector potential power spectrum with the number of 
bins. Again, the change in the number of bins is negligible. 
In this plot, the dashed lines show the power spectra com¬ 
puted with the extra factors of k explicitly included whilst 
summing over the modulus squared values of the held, see 
the velocity consistency check for more information. As ex¬ 
pected, this change affects things the most in the largest bins 
and therefore on large scales and for the smallest number of 
bins, however it does not affect our results. 
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